Subgap current in superconducting tunnel junctions with diffusive electrodes 
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We calculate the subgap current in planar superconducting tunnel junctions with thin-film diffusive leads. It 
is found that the subharmonic gap structure of the tunnel current scales with an effective tunneling transparency 
which may exceed the junction transparency by up to two orders of magnitude depending on the junction ge- 
ometry and the ratio between the coherence length and the elastic scattering length. These results provide an 
alternative explanation of enhanced values of the subgap current in tunneling experiments often ascribed to im- 
perfection of the insulating layer. We also discuss the effect of finite lifetime of quasiparticles as the possible 
origin of additional enhancement of multiparticle tunnel currents. 
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Subgap quasiparticle current in superconducting junctions 
at small applied voltages eV < 2A is the subject of persis- 
tent theoretical interest and experimental research. Recently, 
the problem has attracted new attention, and a number of 
measurements of the subgap current in high-quality tunnel 
junctions have been performed, 1,2 motivated by the problem 
of decoherence in Josephson-junction-based superconducting 
qubits.i The subgap current at zero temperature is due to 
multiparticle tunneling (MPT) processes, 4 whose intensities 
strongly depend on the quality of the insulating layer, being 
enhanced by disorder, localized electronic states, pinholes, 
etci^ The effect of disorder in the junction electrodes on the 
subgap current has never been questioned. 

According to the MPT theory, 4 the subgap tunnel current 
depends on the transparency D of the tunnel barrier: it de- 
creases with decreasing voltage in a steplike fashion with step 
heights proportional to (D/2)" at voltages eV = 2A/n, n = 
1,2... [subharmonic gap structure (SGS)]. Similar results have 
been obtained for junctions with ballistic electrodes^ and 
mesoscopic point contacts with diffusive electrodes 7 on the 
basis of the theory of multiple Andreev reflections (MAR). 5 

Experimentally, the SGS scaling parameter in atomic size 
junctions nicely agrees with the theory; 8 however, in macro- 
scopic tunnel junctions it is usually much larger 1 - 2 (see also 
earlier data 9 ); moreover, there is a smooth residual current 
at a very low voltage. 1 Although enhanced SGS in high- 
transmission junctions could be explained by assuming ran- 
domly distributed resonant levels within the tunnel barrier, 10 
enhanced subgap current in low-transmission junctions with 
presumably good insulating layers remains an open question. 

In this paper we reexamine the problem of the subgap cur- 
rent in macroscopic tunnel junctions, and consider the effects 
of diffusive electrodes and planar junction geometry common 
for the experiment (see Fig.0. Our main result is that the SGS 
scaling parameter for such junctions significantly exceeds the 
junction transparency: for the sandwich-type junction with 
thin-film leads shown in Fig. ^b), the scaling is determined 
by the effective transparency defined as 

D eff ={3^/U)D, (1) 
where E,o = y / D/2A is the diffusive coherence length (we as- 
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FIG. 1: One-dimensional (a) and planar (b) models of the tunnel 
junction with diffusive leads; equilibrium reservoirs are far from the 
contact, L ^> Re- 
sume h = ks = 1), £ is the elastic mean free path, d <C £,o is 
the thickness of the leads, and V is the diffusion coefficient. 
For Al junctions with I ~ d = 50 nm and = 300 nm, the en- 
hancement factor approaches 100. For the junctions with one- 
dimensional (ID) geometry of Fig.^a), D e f = (3^q/£)D. 11 
This result also applies to nonhomogeneous tunnel barriers as 
soon as the size of pinholes (more transparent spots) exceeds 
the elastic mean free path, otherwise the ballistic scaling 6 will 
be valid. 

The enhancement effect can be qualitatively understood 
by considering a short point contact with the reservoirs lo- 
cated very close to the contact, L <C £,o [cf. Fig. [fla)]. In 
this case, the current can be calculated within the mesoscopic 
approach^ 2 " by integrating over contributions of normal con- 
ducting eigenmodes with randomly distributed transparencies. 
The relevant distribution is known to be spread over the in- 
terval ~ [L/£)D ^> D. 13 The most transparent modes domi- 
nate the subgap current, giving D e ff ~ (Ljl)D. In our case 
of junctions with large distance to the reservoirs, the scale of 
the spatial variation of the Green's function <jjo plays the role 
of the effective junction length giving qualitatively our result, 
D e ff ~ (^q/£)D. 14 We note that for the long junctions under 
consideration the statistics of the eigenmode transparencies is 
not known, and a quantitative result has to be derived from the 
quasiclassical theory for diffusive superconductors. 

Our analysis is based on the diffusive equations 15 for the 
quasiclassical two-time Keldysh-Green functions G(r,ti,t2), 

[H, oG] = WW J, GoG = 1, G = ( % V (2) 
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Here H(t\,t2) = [io z d tl + Aexp(/<7,0)/(7 v ]5(fi —tj), (j) is the 
superconducting phase, the sign o denotes time convolution, 
and J = Go VG. Equation Q can be decomposed into the Us- 
adel equation for the retarded or advanced Green's functions 
g R ' and the equation for the Keldysh function G K = g R o / — 
f o g^, where f = f + <7,/_ is the distribution function. 

We present detailed calculations for the simpler, ID geom- 
etry of Fig. ^ a )- At the left electrode, x = —L, the Fourier 
transformations of the two-time functions g and / with respect 
to the variable t\ — ti are given by the equilibrium expressions 

g(E) = a z u(E)+ia y v(E), f(E) = tanh(£/2r), (3) 
(«,v) = (£,A)/§(£), ^ R ' A — [(E ±/0) 2 — A 2 ] 1 / 2 . (4) 

At the right voltage-biased electrode, x = L, the function G is 
defined through the gauge transformatiorU& 

G{L) = G(-L) = S(h)G(-L)S^t 2 ), (5) 

with a unitary operator S(t) = exp[ia z <j) [t)/2], where the phase 
<j> satisfies the Josephson relation 0(f) = 2eVt. 

The boundary conditions 17 for the functions G and J at the 
tunnel barrier (x = ±0) are given by the relations 

/-o = /+o = ^[G-o,oG +0 ], W = ^ = ^|a (6) 

where R is the resistance of the tunnel barrier, Rq — E,o/g is the 
resistance of a piece of the lead with length §o> an d 8 i s me 
conductance of the leads per unit length. Assuming a small 
value of the tunneling parameter W, we neglect the charge 
imbalance function /_ and the superfiuid momentum within 
the leads, as well as the variation of A. In such an approxima- 
tion, Eq. (|5j extends to the whole right lead, G(x) = G(—x) for 
< x < L. The problem is therefore reduced to the solution 
of a static equation for the function G(x,ti,t2) at — L < x < 
with the time-dependent boundary condition (jS) at x = — 0. 
The electric current is related to the Keldysh component J K of 
the matrix current/as 7(f) = (ng /4-e)Tr a z J K (x,t ,f). 15 Using 
Eqs. and 0, it can be expressed as 

I(t) = (n/8eR)TTO z [G,oG) K (t,t). (7) 

In this and following equations, the functions are taken at the 
boundary x = —0. Expanding all functions over harmonics of 
the Josephson frequency, A(E,t) — £ m A(£ ,m)e~ 2leVmt [t = 
(f i + ?2)/2], we arrive at the equation for the dc current /, 

7 = \hk I-J E Tr £"' ^ E ^ {E -"0 

-l(E,m)G K (E,-m)], h = o z f-g^o z . (8) 

In the tunneling limit W <C 1 , the amplitude of the mth har- 
monic is proportional to W m ; thus the zero harmonic m = of 
the functions g and G K in Eq. (|8} plays the key role, while the 
high-order harmonics can be neglected. The effect of these 
harmonics will be discussed later. Within this approximation, 
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FIG. 2: Circuit representation of charge transport in a diffusive tunnel 
junction, eV = 2.5A. 

the Green's function matrix structure in Eq. holds, and the 
current Eq. l|8} exactly transforms to the form 

I=—[ dEN(E)N(E~eV)\n(E-eV)-n(E)}. (9) 
eRJ-oa 

Here N(E) — Re u R is the density of states (DOS) normalized 
to its value in the normal state, and the distribution function 
n = (1/2) (1 — /) approaches the Fermi function hf in equilib- 
rium. Furthermore, we split the integral in Eq. (|9jl into pieces 
of length eV, denoting A k (E) = A(E + keV), 

1 f eV 

/=— J q dEJ(E), J{E)=Y jk= _Jk{E), (10) 
h = ("ft- 1 - n k )pl 1 , p A 7 1 = N k N k - 1 . (11) 

The distribution function n(E,x = —0) satisfies the recurrence 
relation following from the kinetic equation d x [D + d x n) = 0, 

®{\E k \-A)[n F {E k )-n k ] = r(j k+1 -j k ), (12) 

where 0(x) is the Heaviside step function, r = Rn/R-^ 1, and 
Rn is the normal resistance of the lead. To justify Eq. d!2l >. we 
note that the diffusion coefficient D + = (1/2)(1 + \u\ 2 - |v| 2 ) 
is approximately constant, D + w 1, at \E\ > A, which leads 
to the linear function n(E,x) = n_o + (x/L)(n_Q — rip). At 
\E\ < A, D + turns to zero at \x\ 3> §o> which reflects com- 
plete Andreev reflection in the leads and results in zero prob- 
ability current, D+d x ti = 0. Then, using the boundary con- 
dition at the tunnel barrier following from Eq. (|6j, D + d x n = 
(2W/^ )Lk=±iNN k (n k -n), we arrive at Eq. fPH 

A convenient interpretation of Eqs. il It and d!2t in terms 
of circuit theory 18 is given by an infinite network in the en- 
ergy space with the period eV, graphically presented in Fig. [2] 
The electric current spectral density J(E) consists of partial 
currents j k , which flow through the tunnel "resistors" p k con- 
nected to adjacent nodes having "potentials" n k and n k -\. At 
\E\ > A, the nodes are also attached to the distributed "equi- 
librium source" np(E) through equal resistors r. Below we 
impose the equilibrium quasiparticle distribution at \E\ > A, 
n(E) = np(E), neglecting the effect of small resistors r. 

The currents flowing between the nodes outside the gap are 
related to the thermal current; at T = 0, these nodes have equal 
populations (n k = 1 at E k < —A, n k = at E k > A), thus the 
corresponding partial currents are zero, and the thermal cur- 
rent vanishes. As a result, only the current jo flowing across 
the gap through the resistor po survives at T — 0. 

Taking the DOS in the BCS form N = N S = Re(£ we 
see that if any node falls into the gap, the adjacent resistances 
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FIG. 3: DOS and subgap circuits at the applied voltages eV = 1.2A 
(a) and 0.7A (b), for the tunneling parameter W = 10~ 3 . 

turn to infinity, and the current vanishes. For this reason, the 
network period must exceed the gap, eV > 2A, and the integra- 
tion in Eq. d!0> is confined to the region A < E < eV — A. This 
recovers the tunneling model result for the single-particle tun- 
nel current:— the current appears above the threshold, eV = 
2 A, having the threshold value h (2A) = %Aj2eR. 

To evaluate the subgap current, eV < 2A, the DOS must be 
calculated to next order in the parameter W, which requires 
solution of the equations for g following from Eqs. (|2ji and 
(|6j. Using the standard parametrization g = G z e a " 6 ', we obtain 
the equation and the boundary condition for the spectral angle 

e, 

sinh(0-0 s ) = ^081^105, z = x/$q, (13) 
a0 + Wsinh0(cosh0! +cosh0_i) =0 (z = -0). (14) 

With exponential accuracy, the solution of Eq. Jl 3I > for z < 
can be approximated by the formula for a semi-infinite wire, 

tanh{[0( z ) - 9 s ]/4} = tanh[(0_ o - 9 s )/4] exp(fc), (15) 

where k~ l (E) — -y/sinh 9s- Equation dl 51 describes the de- 
cay of perturbations of the spectral functions at distances > <^o 
from the barrier, where the spectral angle approaches its bulk 
value 0$ — arctanh(A/£'). The boundary value of is to be 
found from the equation following from Eqs. i\4\ and 1151 . 

2/fcsinh[(0 5 - 0)/2] = W sinh (cosh 0i + cosh0_i). (16) 

A direct expansion of 9 with respect to W in Eq. d!6i leads to 
the following expression for the DOS within the BCS gap, 

N(E)=W(l -E 2 /A 2 )- 5/4 [N s (E+eV)+N s (E-eV)]. (17) 

The DOS divergencies at \E\ = A, A — eV in Eq. d!7i are po- 
tentially dangerous (cf. Refs. 0), but they can be eliminated 
by improving the perturbation procedure by solving a set of 
recurrences in Eq. (I16> in the vicinity of these points. 

As follows from Eq. (1171 . the tunneling processes transfer 
the DOS in the energy space into the BCS gap at the distances 
±eV from the regions \E\ > A, thus forming an effective spa- 
tial potential well of the width ~ <^o at the tunnel barrier. At 



eV > A the BCS gap is entirely filled with the quasiparticle 
states with a small local DOS ~ W, as shown in Fig. 0a). 
The appearance of localized states enables the quasiparticles 
to overcome the BCS gap at eV < 2A via two steps involv- 
ing intermediate Andreev reflection at energies \E\ < A. The 
population of the intermediate state cannot be taken to be in 
equilibrium because the subgap quasiparticles cannot access 
the equilibrium electrodes. In the circuit terms, the node k = 
is disconnected from the equilibrium source, and the subgap 
current flows through the two large resistances po,Pi ~ W~ l 
(two-particle current), see Fig.QJa). The corresponding partial 
currents are equal, 7*0 = j\ = [np(E\) — nf (i?-i)]/(po + Pi), 
and their contribution to I(V) is confined to the energy region 
< E < eV — A (a similar contribution at A < E < eV comes 
from jo and j-i). Thus the two-particle current appears above 
the threshold eV — A, having the threshold value /2(A) = 
TtWA/eR = 2W/j(2A). At eV = 2 A, the two-particle current 
exhibits a sharp peak with the height /2(2A) » 23W 2 l 5 A/eR\ 
at larger voltages, it approaches a constant value giving rise to 
the excess current I exc w 6.2W 2 ^A/eR. 

At eV < A, a minigap opens in the DOS around the zero 
energy [see Fig.]3jb)], however, since the number of subgap 
resistors increases up to three (three-particle current), the cur- 
rent across the minigap will persist as long as the network 
period exceeds the minigap size, eV > 2(A — eV), i.e., at 
eV > 2A/3. The central resistance po is large, po ~ W~ 2 , and 
dominates the net subgap resistance. This leads to a smaller 
charge current with the threshold value ^(2 A/3) ~ 2W/2(A). 
At eV < 2A/3 the network period becomes smaller than the 
minigap, and further correction to the DOS is required. 

Similar results were found for the planar junction Fig.^b), 
using the equation for the functions G±q at the top (+0) and 
bottom (—0) sides of the tunnel barrier, i\a z E + ia y A, G_o] = 
2AW[G_o,G ; +o], with the modified tunneling parameter W = 
(3^q /4£d)D. This equation is derived by averaging Eq. Q 
over the thickness of overlapping leads and using Eq. l|6} (cf. 
Ref.Eol). From this equation we obtain a relation for the spec- 
tral angle that does not significantly differ from Eq. dl 6i . 

£ 2 sinh(0,s- 0) = Wsinh0(cosh0i +cosh0_i), (18) 

thus giving results which are close to those for the ID model 
with the same magnitude of the parameter W. 

The presented calculation scheme, combining circuit theory 
arguments with DOS iteration procedures, suggests an appeal- 
ingly simple explanation for the diffusive SGS: the decreasing 
applied voltage results in a shrinking period of the network in 
Fig. |2 hence a stepwise increase of the number of subgap re- 
sistors involved; simultaneously, the number of DOS steps, 
scaled as W, increases, as shown in Fig. Ea). This results 
in the current staircase with the height of the steps given by 
/„ - (2W)"- 1 /!, at 2A/n < eV < 2A/(n - 1). The quantitative 
result for the current at arbitrary voltages and temperatures is 



I(V) = 



eV 



dE N+ +N- . . [°° 2dE , 

{n--n+)+ — — (n F -n F i), 
J A eRpi 



/o eR Pa J a eRp\ 

N ± =lnt[{ATE)/eV] + l, n ± (E) = n F (E ±N± ). (19) 

In this equation, the second term represents the thermal cur- 
rent, the integers ±N± are the indices of the nodes closest to 
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FIG. 4: DOS at eV = 0.4A (a) and I-V characteristics (b) for the 
tunneling parameter W = 10~ 3 and two values of the damping pa- 
rameter: 7 = (solid line) and y = 0.003A (dashed line). 

the gap edges outside the gap, Int(x) denotes integer part of x, 
and the quantity Pa(E) = y* k =\_ N Pk has the meaning of net 
subgap resistance. The subgap distribution function reads 

»(£)=« + + («_-» + )p A 1 £^ 1 p,. (20) 

Equations J19I and \20\ are the main technical results of 
the paper. The I-V characteristic (IVC) of the planar tun- 
nel junction calculated from Eqs. dl 91 and Jl 8i at T = and 
shown in Fig. |4jb), was found to be very close to the result 
for a ballistic point contact 6 with the effective transparency 
D e ff — 4W — (3^q /ld)D. This justifies our statement made in 
the introduction, and is the main conclusion of this paper. 

In low-transmissive junctions, enhanced subgap current at 
eV < A has been observed (see, e.g., Ref.Q])- This anomaly 
might be due to many-body interaction effects which intro- 
duce a finite lifetime (damping) of the quasiparticles. The 
damping effect can be qualitatively modeled by a small imag- 
inary addition to the energy in the spectral functions, E — > 
E + iy. This would lead to a small residual DOS within the 
BCS gap and cut the DOS staircase at the level of the order 
of y/A, see Fig. |4ja). This will result in the smearing of the 



tunneling SGS and crossover to a linear IVC at low voltages, 
/ = 2.2(y/A) 2 V /R, similar to the incoherent MAR regime. 18 
The IVC calculated from Eq. {19} for y = 0.003A and shown 
in Fig.|4lb) by a dashed line confirms these considerations. 

We conclude our analysis with the estimation of the contri- 
bution of higher harmonics of the functions g and G K to the dc 
current. At T = 0, the contribution Si of the first harmonics 
\m\ = 1 (the higher harmonics, \m\ > 1, are smaller, ~ W m ) is 



81= / c/£lmvlm( -cosh 2 - + -sinh 2 - J, (21) 

eR Jo \p 2 q 2 J 

where % = ft + fl_ lt % = Bi + B^, (p,q) 2 = (§f + %*f)/2iA, 
and v = sinh 9. At eV < A, the energy £_ i appears in the sub- 
gap region, where 0*j = 0_i + %i and £, A l = £, R i ; for this 
reason, 81 turns to zero at eV < A, similar to h- Numerical 
calculations show that the contribution of the first harmonics 
to the IVC does not exceed 30%. From this we conclude that 
the adopted quasistatic approach gives a rather good approxi- 
mation to a complete solution. 

In our treatment, we have neglected inelastic scattering, 
which might affect the quasiparticle distribution at subgap en- 
ergies. Analysis shows that this effect becomes essential un- 
der the condition Wt £ A< 1, where T £ is the relaxation time. 
However, this does not affect the estimate of the effective scal- 
ing factor and only changes the details of the IVC shape. 

In conclusion, we have developed a theory of subgap charge 
transport and subharmonic gap structure in superconducting 
tunnel junctions with planar geometry and diffusive thin-film 
electrodes. We found that the role of scaling factor in the sub- 
harmonic gap structure is played by the effective tunneling 
transparency D e g = (3^q / £d)D, which may greatly exceed the 
bare transparency D of the junction tunnel barrier. 
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